{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "c19e838a-0ed1-47ff-a397-ea0053ee6da7", "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "%matplotlib inline\n", "%config InlineBackend.figure_format = 'svg'\n", "\n", "import warnings \n", "warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n", "\n", "import numpy as np\n", "\n", "from mpi4py import MPI\n", "\n", "from h5 import HDFArchive\n", "from itertools import product\n", "\n", "from triqs.plot.mpl_interface import plt, oplot, oplotr, oploti\n", "plt.rcParams[\"figure.figsize\"] = (6,4) # set default size for all figures" ] }, { "cell_type": "markdown", "id": "9738461c-04a1-443e-ae57-8e88d4c193f8", "metadata": {}, "source": [ "# Metal insulator transition\n", "\n", "The half-filled Hubbard model undergoes a transition from a metallic state to an (Mott) insulating state when increasing the interaction strength $U$. One model where this can be solved exactly is the Hubbard model on the Bethe graph in the limit of infinite connectivity. In this limit Dynamical Mean Field Theory (DMFT) is exact and can be used to study the transition.\n", "\n", "The DMFT lattice self-consistency relation on the Bethe graph has the form\n", "\n", "$$\n", "\\Delta_\\sigma(\\tau) = t^2 G_\\sigma(\\tau)\n", "$$\n", "\n", "where $\\Delta_\\sigma(\\tau)$ is the hybridization function, $G_\\sigma(\\tau)$ is the local single particle Green's function with spin $\\sigma \\in \\{\\uparrow, \\downarrow\\}$, and we set the nearest neighbour hoppint $t$ to unity ($t=1$).\n", "\n", "At low temperature the transition is first order, but at higher temperatures it turns into a cross over. Here we will study this crossover at inverse temperature $\\beta = 5$." ] }, { "cell_type": "markdown", "id": "6aa5e06a-d60a-4184-ae2d-b820b0cbf519", "metadata": {}, "source": [ "## Solver setup\n", "\n", "Since we will be performing self-consistent DMFT calculations we modularize the solver setup using a helper function `setup_solver()`." ] }, { "cell_type": "code", "execution_count": 2, "id": "9372debf-a505-423b-a0f6-5d517fc0b3c7", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Warning: could not identify MPI environment!\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Starting serial run at: 2025-08-19 17:41:44.159873\n" ] } ], "source": [ "from triqs_soehyb import Solver\n", "\n", "from triqs.gf import make_gf_dlr_imtime, make_gf_dlr_imfreq, SemiCircular\n", "\n", "def setup_solver():\n", " \n", " S = Solver(\n", " beta=5.0, # Inverse temperature\n", " gf_struct=[('up', 1), ('do', 1)], # Green's function structure 1st index: name, 2nd index: dimension of subspace\n", " eps=1e-12, # Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions\n", " w_max=10.0, # DLR freqiency cut-off (the spectrum of the model must be in the range [-w_max, +w_max]\n", " )\n", " \n", " for bidx, delta_tau in S.Delta_tau:\n", " delta_w = make_gf_dlr_imfreq(delta_tau)\n", " delta_w << SemiCircular(2.0)\n", " delta_tau[:] = make_gf_dlr_imtime(delta_w)\n", "\n", " return S" ] }, { "cell_type": "markdown", "id": "0e76af5c-55eb-48c1-a4d7-dee2b28a4793", "metadata": {}, "source": [ "## Local Hubbard interaction\n", "\n", "The local interaction in the Hubbard model has the form\n", "\n", "$$\n", "H_{int} = U n_\\uparrow n_\\downarrow - \\mu ( n_\\uparrow + n_\\downarrow )\n", "$$\n", "\n", "where $U$ is the so called Hubbard-$U$ and $\\mu$ is the chemical potential, which at half filling is given by $\\mu = U/2$. The helper function `get_h_int(U)` constructs the local Hamiltonian using Triqs second-quantization operators." ] }, { "cell_type": "code", "execution_count": 3, "id": "4359cae7-6969-4ad4-9b12-67b5eb707d70", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "h_int = -2.5*c_dag('do',0)*c('do',0) + -2.5*c_dag('up',0)*c('up',0) + 5*c_dag('do',0)*c_dag('up',0)*c('up',0)*c('do',0)\n" ] } ], "source": [ "from triqs.operators import n\n", "\n", "def get_h_int(U, mu=None):\n", " if mu is None:\n", " mu = U / 2\n", " h_int = U * n('up',0) * n('do', 0) - mu * ( n('up', 0) + n('do', 0) )\n", " return h_int\n", "\n", "h_int = get_h_int(U=5.0)\n", "\n", "print(f'h_int = {h_int}')" ] }, { "cell_type": "markdown", "id": "71edb549-fc4e-40a0-ba27-c6da437a6989", "metadata": {}, "source": [ "## DMFT self-consistency\n", "\n", "Finally we design a function `solve_dmft_self_consistency(S, h_int, order, ...)` that takes a solver `S`, a local Hamiltonian `h_int`, the wanted expansion `order`, and performs the DMFT self-consistency iterations." ] }, { "cell_type": "code", "execution_count": 4, "id": "0f9aa03f-15a5-4480-a5e2-ce0d99290d9d", "metadata": {}, "outputs": [], "source": [ "def solve_dmft_self_consistency(S, h_int, order, \n", " dmft_maxiter=100, dmft_tol=1e-6):\n", "\n", " for dmft_iter in range(1, dmft_maxiter+1):\n", "\n", " S.solve(h_int=h_int, order=order, maxiter=100, verbose=False)\n", " \n", " diff = np.max(np.abs(S.Delta_tau['up'].data - S.G_tau['up'].data))\n", "\n", " print(f'DMFT: dmft_iter = {dmft_iter}, diff = {diff:2.2E}')\n", " \n", " if diff < dmft_tol: break\n", " \n", " S.Delta_tau[:] = S.G_tau # DMFT Bethe lattice self consistency\n", "\n", " S.d = S.S.get_expectation_value(n('up', 0) * n('do', 0)).real\n", " S.n_up = S.S.get_expectation_value(n('up', 0)).real\n", " S.n_do = S.S.get_expectation_value(n('do', 0)).real\n", "\n", " return S" ] }, { "cell_type": "markdown", "id": "6dc1313a-0310-492c-a9b0-02555ac952b7", "metadata": {}, "source": [ "## Solve NCA and OCA for a sweep in $U$\n", "\n", "To study the metal-insulator cross over we solve the Hubbard model on the Bethe graph for a range of $U$ values at first and second expansion order (a.k.a. the non-crossing and one-crossing approximation NCA/OCA)." ] }, { "cell_type": "code", "execution_count": 5, "id": "66c1de8d-a6e5-4eb1-9d95-ed8b0fd41d70", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " ___ ___ _ ___ _____\n", " / __| ___| __|__| || \\ \\ / / _ )\n", " \\__ \\/ _ \\ _|___| __ |\\ V /| _ \\\n", " |___/\\___/___| |_||_| |_| |___/ [github.com/TRIQS/soehyb]\n", "\n", "beta = 5.0, lamb = 5.00E+01, eps = 1.00E-12, N_DLR = 27\n", "fundamental_operators = [1*c('up',0), 1*c('do',0)]\n", "H_mat.shape = (4, 4)\n", "H_loc = 0\n", "\n", "--> order = 1, U = 3.0\n", "DMFT: dmft_iter = 1, diff = 1.27E-01\n", "DMFT: dmft_iter = 2, diff = 1.24E-02\n", "DMFT: dmft_iter = 3, diff = 1.58E-03\n", "DMFT: dmft_iter = 4, diff = 3.43E-04\n", "DMFT: dmft_iter = 5, diff = 6.34E-05\n", "DMFT: dmft_iter = 6, diff = 1.24E-05\n", "DMFT: dmft_iter = 7, diff = 2.39E-06\n", "DMFT: dmft_iter = 8, diff = 4.62E-07\n", "--> Storing: data_order_1_U_03.00000000.h5\n", "--> order = 1, U = 4.0\n", "DMFT: dmft_iter = 1, diff = 4.31E-02\n", "DMFT: dmft_iter = 2, diff = 1.11E-02\n", "DMFT: dmft_iter = 3, diff = 3.24E-03\n", "DMFT: dmft_iter = 4, diff = 9.52E-04\n", "DMFT: dmft_iter = 5, diff = 2.83E-04\n", "DMFT: dmft_iter = 6, diff = 8.44E-05\n", "DMFT: dmft_iter = 7, diff = 2.52E-05\n", "DMFT: dmft_iter = 8, diff = 7.50E-06\n", "DMFT: dmft_iter = 9, diff = 2.24E-06\n", "DMFT: dmft_iter = 10, diff = 6.68E-07\n", "--> Storing: data_order_1_U_04.00000000.h5\n", "--> order = 1, U = 5.0\n", "DMFT: dmft_iter = 1, diff = 3.63E-02\n", "DMFT: dmft_iter = 2, diff = 5.84E-03\n", "DMFT: dmft_iter = 3, diff = 1.63E-03\n", "DMFT: dmft_iter = 4, diff = 4.16E-04\n", "DMFT: dmft_iter = 5, diff = 1.06E-04\n", "DMFT: dmft_iter = 6, diff = 2.69E-05\n", "DMFT: dmft_iter = 7, diff = 6.85E-06\n", "DMFT: dmft_iter = 8, diff = 1.74E-06\n", "DMFT: dmft_iter = 9, diff = 4.42E-07\n", "--> Storing: data_order_1_U_05.00000000.h5\n", "--> order = 1, U = 6.0\n", "DMFT: dmft_iter = 1, diff = 3.04E-02\n", "DMFT: dmft_iter = 2, diff = 2.63E-03\n", "DMFT: dmft_iter = 3, diff = 5.08E-04\n", "DMFT: dmft_iter = 4, diff = 1.05E-04\n", "DMFT: dmft_iter = 5, diff = 2.09E-05\n", "DMFT: dmft_iter = 6, diff = 4.09E-06\n", "DMFT: dmft_iter = 7, diff = 8.00E-07\n", "--> Storing: data_order_1_U_06.00000000.h5\n", " ___ ___ _ ___ _____\n", " / __| ___| __|__| || \\ \\ / / _ )\n", " \\__ \\/ _ \\ _|___| __ |\\ V /| _ \\\n", " |___/\\___/___| |_||_| |_| |___/ [github.com/TRIQS/soehyb]\n", "\n", "beta = 5.0, lamb = 5.00E+01, eps = 1.00E-12, N_DLR = 27\n", "fundamental_operators = [1*c('up',0), 1*c('do',0)]\n", "H_mat.shape = (4, 4)\n", "H_loc = 0\n", "\n", "--> order = 2, U = 3.0\n", "DMFT: dmft_iter = 1, diff = 8.02E-02\n", "DMFT: dmft_iter = 2, diff = 3.40E-03\n", "DMFT: dmft_iter = 3, diff = 3.65E-04\n", "DMFT: dmft_iter = 4, diff = 5.09E-05\n", "DMFT: dmft_iter = 5, diff = 7.03E-06\n", "DMFT: dmft_iter = 6, diff = 9.69E-07\n", "--> Storing: data_order_2_U_03.00000000.h5\n", "--> order = 2, U = 4.0\n", "DMFT: dmft_iter = 1, diff = 4.58E-02\n", "DMFT: dmft_iter = 2, diff = 1.46E-02\n", "DMFT: dmft_iter = 3, diff = 6.23E-03\n", "DMFT: dmft_iter = 4, diff = 2.71E-03\n", "DMFT: dmft_iter = 5, diff = 1.20E-03\n", "DMFT: dmft_iter = 6, diff = 5.39E-04\n", "DMFT: dmft_iter = 7, diff = 2.42E-04\n", "DMFT: dmft_iter = 8, diff = 1.09E-04\n", "DMFT: dmft_iter = 9, diff = 4.91E-05\n", "DMFT: dmft_iter = 10, diff = 2.22E-05\n", "DMFT: dmft_iter = 11, diff = 9.99E-06\n", "DMFT: dmft_iter = 12, diff = 4.50E-06\n", "DMFT: dmft_iter = 13, diff = 2.03E-06\n", "DMFT: dmft_iter = 14, diff = 9.16E-07\n", "--> Storing: data_order_2_U_04.00000000.h5\n", "--> order = 2, U = 5.0\n", "DMFT: dmft_iter = 1, diff = 4.13E-02\n", "DMFT: dmft_iter = 2, diff = 1.34E-02\n", "DMFT: dmft_iter = 3, diff = 6.27E-03\n", "DMFT: dmft_iter = 4, diff = 2.80E-03\n", "DMFT: dmft_iter = 5, diff = 1.23E-03\n", "DMFT: dmft_iter = 6, diff = 5.41E-04\n", "DMFT: dmft_iter = 7, diff = 2.37E-04\n", "DMFT: dmft_iter = 8, diff = 1.04E-04\n", "DMFT: dmft_iter = 9, diff = 4.56E-05\n", "DMFT: dmft_iter = 10, diff = 2.00E-05\n", "DMFT: dmft_iter = 11, diff = 8.75E-06\n", "DMFT: dmft_iter = 12, diff = 3.84E-06\n", "DMFT: dmft_iter = 13, diff = 1.68E-06\n", "DMFT: dmft_iter = 14, diff = 7.37E-07\n", "--> Storing: data_order_2_U_05.00000000.h5\n", "--> order = 2, U = 6.0\n", "DMFT: dmft_iter = 1, diff = 3.30E-02\n", "DMFT: dmft_iter = 2, diff = 5.61E-03\n", "DMFT: dmft_iter = 3, diff = 1.83E-03\n", "DMFT: dmft_iter = 4, diff = 6.47E-04\n", "DMFT: dmft_iter = 5, diff = 2.18E-04\n", "DMFT: dmft_iter = 6, diff = 7.20E-05\n", "DMFT: dmft_iter = 7, diff = 2.37E-05\n", "DMFT: dmft_iter = 8, diff = 7.77E-06\n", "DMFT: dmft_iter = 9, diff = 2.55E-06\n", "DMFT: dmft_iter = 10, diff = 8.34E-07\n", "--> Storing: data_order_2_U_06.00000000.h5\n" ] } ], "source": [ "orders = [1, 2]\n", "Us = np.arange(3.0, 7.0, 1.0)\n", "\n", "for order in orders:\n", "\n", " S = setup_solver()\n", " \n", " for U in Us:\n", " \n", " print(f'--> order = {order}, U = {U}')\n", " \n", " S.U = U\n", " h_int = get_h_int(U)\n", "\n", " S = solve_dmft_self_consistency(S, h_int, order)\n", "\n", " filename = f'data_order_{order}_U_{U:011.8f}.h5'\n", " print(f'--> Storing: {filename}')\n", " with HDFArchive(filename, 'w') as A: \n", " A['S'] = S" ] }, { "cell_type": "markdown", "id": "12329eef-af0d-46b1-a48e-98a10798ed70", "metadata": {}, "source": [ "## Load result\n", "\n", "To visualize the crossover we load the results from file using a helper function and class." ] }, { "cell_type": "code", "execution_count": 6, "id": "4b274fc5-df67-4570-8a2f-a2dad74c7fc3", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--> Loading: data_order_1_U_03.00000000.h5\n", "--> Loading: data_order_1_U_04.00000000.h5\n", "--> Loading: data_order_1_U_05.00000000.h5\n", "--> Loading: data_order_1_U_06.00000000.h5\n", "--> Loading: data_order_2_U_03.00000000.h5\n", "--> Loading: data_order_2_U_04.00000000.h5\n", "--> Loading: data_order_2_U_05.00000000.h5\n", "--> Loading: data_order_2_U_06.00000000.h5\n" ] } ], "source": [ "import glob\n", "\n", "class AttrList:\n", " def __init__(self, l): self.l = l\n", " def __getattr__(self, attr): return [ getattr(x, attr) for x in self.l ]\n", "\n", "\n", "def load_order(order):\n", "\n", " filenames = np.sort(glob.glob(f'data_order_{order}_U_*.h5'))\n", "\n", " results = []\n", " for filename in filenames:\n", " print(f'--> Loading: {filename}')\n", " with HDFArchive(filename, 'r') as A:\n", " results.append(A['S'])\n", "\n", " return AttrList(results)\n", "\n", "nca = load_order(1)\n", "oca = load_order(2)" ] }, { "cell_type": "markdown", "id": "95e8bf34-798a-47de-8186-b41125588626", "metadata": {}, "source": [ "## Visualization\n", "\n", "To visualize the metal insulator cross over we plot the double occupancy $\\langle n_\\uparrow n_\\downarrow \\rangle$ as a function of $U$ for NCA and OCA." ] }, { "cell_type": "code", "execution_count": 7, "id": "a8d689c5-71b3-4fed-be16-937572394579", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2025-08-19T17:42:11.588108\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.3, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(nca.U, nca.d, 'o-', label='NCA')\n", "plt.plot(oca.U, oca.d, 'o-', label='OCA')\n", "\n", "plt.ylabel(r'$\\langle n_{\\uparrow} n_{\\downarrow} \\rangle$')\n", "plt.xlabel('U')\n", "plt.ylim(bottom=0, top=0.12)\n", "plt.legend(loc='best');" ] }, { "cell_type": "markdown", "id": "91b02ba5-6e6d-488a-90ff-22e96f9ccea1", "metadata": {}, "source": [ "This reproduces the results in Fig. 4a in [M. Eckstein, P. Werner, Phys. Rev. B 82, 115115 (2010)](https://doi.org/10.1103/PhysRevB.82.115115)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.17" } }, "nbformat": 4, "nbformat_minor": 5 }